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Abstract 

I calculate the hydrodynamic limit of the BGK approximation of the 
Boltzmann equation for the case of a long stress relaxation time and find 
that the stress obeys a viscoelastic constitutive equation. The constitu- 
tive equation is different from standard constitutive equations used for 
polymeric liquids in that it is not "objective" because inertial effects are 
important. I calculate the exact solution for the stresses for simple shear 
flow and elongational flow. 

1 Introduction 

It has long been known that not only complex fluids like polymeric liquids show 
viscoelastic behavior but that even gases show some viscoelasticity Q . Because 
these effects are small for gases and difficult to measure in most cases they 
have only received a more detailed consideration when molecular dynamics sim- 
ulations showed non-Newtonian behavior in simple shear flow and Zwanzig 
H found an analytical solution for the BGK approximation of the Boltzmann 
equation that reproduced the shear thinning found in the computer simulations. 
Examination of the non-Newtonian aspects of flows described by the Boltzmann 
equation focuses on the case of a simple shear flow [|[ |], |(| 0] and it was analyzed 
in increasing detail by generalizing the results to hard spheres [|[ [|. Recently 
the results have been generalized to a perturbative expansion around the simple 
shear flow by Lee et al. 0. 

Giraud and d'Humieres|^| take a different angle on the non-Newtonian limit 
of the Boltzmann equation. They have attempted to use the viscoelastic proper- 
ties hidden in the BGK-Boltzmann equation to define lattice Boltzmann models 
with viscoelastic properties. So far, however, their results are limited to linear 
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viscoelasticity. We hope that the analytical methods derived in this article can 
help to direct the research into viscoelastic lattice Boltzmann methods. 

In this article we present an intriguingly simple analysis of the moments of 
the Boltzmann equation in the BGK approximation which allows us to derive 
a viscoelastic constitutive equation. This constitutive equation has a number 
of interesting properties. Firstly it is not "objective", i.e. inertial effects can 
not be neglected. Secondly the first normal stress difference has a negative 
sign as compared with a positive sign for polymeric liquids. Unlike the upper- 
convected Maxwell model used to describe polymeric liquids our constitutive 
equation shows shear thinning. 

In the second part of this article we will show that our constitutive equation 
reduces to the analytical results derived by Zwanzig[|| for the case of simple 
shear flow. For this simple situation we then give an intuitive explanation for 
the difference in viscoelastic behavior for a gas and a polymeric liquid. We also 
find a new exact solution for the stress for elongational flows. 

2 Derivation of viscoelasticity in the BGK Boltz- 
mann equation 

The Boltzmann equation was derived for rarefied gases and is a complicated 
integral equation. However in the hydrodynamic limit the Navier-Stokes equa- 
tions for fluid flow can be recovered and this is one of the reasons that it is 
also the basis for a popular method of simulating fluids, the lattice Boltzmann 
method. All the properties required for the description of a fluid are retained 
in the somewhat simpler single relaxation time approximation introduced by 
Bhatnagar et al. (BGK) which is still a complicated non-linear equation. In 
certain simple situations, however, this equation has an analytical solution jfl]. 

In this section we derive a viscoelastic limit of the BGK approximation Hof 
the Boltzmann equation in the limit of large relaxation times. The Boltzmann 
equation is an evolution equation for probability density f v (x, t) to find a par- 
ticle at time t, at position x with velocity v. The evolution equation consists of 
a free streaming of the particles and a collision term. The effect of the collision 
term is approximated in the BGK approximation as a simple relaxation of the 
distribution function towards the equilibrium distribution / . This local equi- 
librium distribution is the Maxwell-Boltzmann distribution for a given density, 
velocity and temperature. The evolution equation is 

d t fv + v a d a f v = -(f°-f v ) (1) 

T 

where the relaxation time r is r(p,9) = To/(pT b ). The exponent b can be 
calculated for molecules with potential V(r) ~ r~ l to be b — 1/2 — 2/1 ||. The 
simplest case is the case of so called Maxwell molecules (I = 4) for which 6 = 
and there is not T-dependence. We define the density p, the mean velocity u a 
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and the temperature 9 as 

J f v dv = p, J f v v a = pu a , J /„(v - u) 2 = Dp8 (2) 

where D is the number of spatial dimensions. Because of mass, momentum and 
energy conservation the equilibrium distribution has to have the same moments 

f°dv = p, f"v a =pu a , f"(v a - u a )(vp - up) = p65 a p (3) 



and the second moment is isotropic. We can now integrate the Boltzmann equa- 
tion over the conserved moments 1, v, v 2 corresponding to mass, momentum and 
energy and get the continuity equation 

d t p + d a {pu a ) =0 (4) 

and a momentum conservation equation 

pd t u a + pupdpu a = -dpTLap (5) 

where 

IIa/3 = / fv(Va ~ U a )(Vj3 - Up)dv (6) 



We now need to calculate Tl a p and we split the contribution to Tl a p into an 
equilibrium and a traceless non-equilibrium part where we use 

f v = f°-T(d t f v +V a d a f v ) (7) 

from equation (|l|) . 

U a p = f f v (v a - u a )(vp - up)dv 

[fv - T (dtfv + d 7 v 7 f v )] (v a - u a )(vp - up) 

= p98 a p + a a p (8) 

where a a p is a traceless stress. 

The energy conservation equation is then 



= (d t f M +v dpf M )dv 

= d t (U aa + pu a u a ) 

+dp(Qp aa + Tlp a u a + U a pU a + Tl aa Up + pupu a u a ) (9) 

where Q is related to the third velocity moment of the distribution function 



Q a p~ ( = I f v (v a - u a )(vp - up)(v-f - u~,)dv (10) 
f v V a VpVjdv — u a Ilp 7 — UpH al — u 7 H a p — pu a upu-y. 
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We can restate the energy conservation equation (g) as 

2 + D 1 2 

d t (p9) +updp(p9) = j^—p9dpup-—dpQp aa -—dpu a a a p (11) 

The equation (|J) imposes a constitutive equation for the stress a a p'- 

-^°af} = ~ J (dtfv + V 7 dyf v )(v a - U a )(V/3 - Up) 

= -d t (p96 a p + a a p) - d 1 up(p95 ai + a ai ) 

-d 1 U a (p95 1 p + (T 7/ a) - 8.yU 7 (p9d a + (Tap) 

-u 1 d 1 {p95 a p + a a p) - d 1 Q a p 1 (12) 
Using the equation for the energy conservation we get for the 9 terms 

dt(p0)5 a p + u 1 d 1 (p6)5 a p + pd(dpu a + d a up + d^u^Sap) 

2 2 1 

= p9(dpu a + d a up - —d-fU^Sap) + — dyUgo-rys - j^d^Q^ss- (13) 

Introducing a convected time derivative as 

&ap<i> = D t (T a p + (d 7 u a a y p + a ai d~fUp), (14) 

where the total derivative is D t = dt + v a d a , we can now write the constitutive 
equation as 

2 

Ca/3 + tct q/ 3<i> - r —(crysdsu^Sap + Td y u 7 a a p 

2 1 

= -Tp8(d a up + dpu a - —d 7 u^S a p) - rd^Qapj + r—d y Q lS sKp (15) 

This expression is similar to the stress relaxation equation proposed by Maxwell [Q 
and we get his result if we replace the convected derivative with a partial time 
derivative. Eqn. ( p"5| ) was first derived in its complete form from an expansion 
of the full Boltzmann equation in Hermite polynomials by Grad[^0[. 

The usual expansion of the Boltzmann equation in the hydrodynamic limit 
is done under the assumption r = 0(1) and d = O(e). Then we get for the 
leading order of the stress (of order e) 

a = -t P 9 (Vu + (Vu) 1 - ^tr(Vu)S^ (16) 

For larger relaxation times, however, we obtain a viscoelastic result. We keep the 
assumption that derivatives are small to order epsilon d = 0(e) but keep terms 
up to order 0(e 2 ). We now get to leading order 0(e 2 ) the full equation ( p"5| ) 
but the terms containing Q still need to be expressed in terms of macroscopic 
quantities. The contribution of the equilibrium distribution to Q vanishes so it 
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is useful to express Q in terms of the first order non-equilibrium contributions. 
If we iteratively substitute eqn. (Q) twice into eqn.([lO|) we can write 

Qapj = -T yj d t fy(v a - U a )(Vp - Up)(v y - 1i 7 ) 

+ J dsfiVs(v a - u a )[vp - m/j)(w 7 - u 7 )^j + 0(d 2 ) 

ffd t [(v a - u a )(vp - up)(v-f - u 7 )] 
+ / fiVsds [{v a - u a )(vp - up)(v~, - it 7 )] 

-ds J f°(v a - u a )(vp - up)(v 7 - U-y)(vs - ug)\ + 0(d 2 ) 

= —t {{5 a p8 1 s + SajSps + 5 a sSp 1 )ds(p9 2 ) 

+ (d t u a + usd8U a )p05pj + {d t up + usdsup)p95 aj 
+(d t ti 7 + u s d s u^)p98 ol p} + 0(d 2 ) 

= -t {(SapS^s + Sa-ySps + 5 a s5p 1 )ds(p9 2 ) 

~{S a pS lS + d aj 5p S + S aS 6p y )6d s (p0)} + 0{d 2 ) 

= -Tpe(s a pS lS + s ai Sps + s Q sSp 7 )d s e + o{d 2 ) (17) 

which is the only approximation we have to make. We then get for the Q terms 
in eqn. ( |l6|) 

= -d a (rp9dpd) ~ dp(T P 9d a 6) + ^d 1 {T P 9d 1 9)5 a p + 0(e 2 ) 

= -rp9^2d a dp~^d 1 d 1 S a p^+0(e 2 ) (18) 

where we have used r(p0) — r /(p9) in the last line. So the viscoelastic consti- 
tutive equation is 

2 

o a p + Ta a p K1> = rp9(d a up + dpu a - — <9 7 "u 7 (5 Q/3 ) 



-rp9 ( 2d a dp - -^9 7 a 7 ) 9 + 0(e 3 ) ( I <) i 



where the last term in 9 represents stresses induced by second derivatives of 
the temperature. Note that this constitutive equation is exact if second order 
derivatives vanish. 
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3 Comparison to usual "objective" constitutive 
equations 

In the theory of viscoelastic constitutive equations usually the assumption is 
made that the constitutive equation should be invariant under an arbitrary 
transformation that preserves distances and time intervals [fll| . That also 
requires form-invariance of the constitutive equation under some non-inertial 
transformations such as rotations and accelerated systems. The underlying idea 
is that the stresses in a polymeric fluid only the stretching of the system is 
relevant and all inertial contributions can be neglected. In particular a con- 
verted time derivative that obeys this symmetry is said to be an "objective 
time derivative" . 

There are two such derivatives; the upper convected derivative 

°"q/3(1) = D t <J al3 — (9 7 M Q tT 7/ 9 + (TaydjUp) (20) 

and the lower convected derivative 

Cal = D tcr a p + {d a u 1 a 1 p + cr a ^dpu^). (21) 

Any linear combination of these two is also objective, and in the past the Jau- 
mann derivative given by (c a| 8(i) + er^)/2 has been popular. Today, however, 
constitutive equations for viscoelastic fluids are constructed with a broad pref- 
erence for the upper convected derivative because these agree best with the new 
experimental results jll| . 

Many frequently used constitutive models are contained in the rather general 
Oldroyd 8-constant model |ll], 7.3-2]. If we only keep constants that can be 
related to our constitutive equation we get 

a aj3 + \icr a f}(i) + -Md-fUgasydap = -r} (d a U/3 + dpu a ) (22) 

as a special case of the Oldroyd 8-constant model. This is identical to our 
constitutive model of eqn. (H) up to the sign change between ct q/ 3(i) and cr a /3<i> 
if we consider that the Oldroyd model is only derived for incompressible fluids 
(<9 7 it 7 = 0) and does not consider elastic effects of temperature gradients (86 — 
0). It is interesting to note that the most frequently used special cases of the 
Oldroyd 8-constant model all assume that \e — 0. (The equation above with 
A6 = is known as the upper convected Maxwell model). 

The difference between the Oldroyd model and our constitutive equation lies 
purely in the nature of the convected derivatives of eqn. ( pp| ) and er<i> of 
eqn. ( |l4| ) which is mathematically a difference in the sign of the dua terms. In 
particular this means that the new convected derivative is not "objective" . It 
is not too surprising, however, that inertial effects are important for a gas and 
therefore there is no reason to require an objective derivative for the viscoelastic 
constitutive equation of a gas. We will be discussing the effect of this subtle 
difference for two special cases, the simple shear flow and the elongational flow 
and give also an intuitive explanation for the differences in the next section. 
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3.1 Simple shear flow 

The most studied situation is probably the simple shear flow with velocity profile 

iv 

1 (23) 


where 7 is the shear rate. The stress is then given by 

@xx &xy 

a X y <Jyy ) (24) 

a zz 

and we obtain the differential equations for the stress terms 

<Jxx + Td t (Txx = t Q - 2^ ja xy (25) 
2 

u vv + rdtOyy = r-ja xy (26) 

a xy + Td t a xy + Tjo-yy = -rpO^ (27) 

To close this system of equations we also need the heat equation (|ll|) which 
reads for a spatially constant temperature and a divergence free velocity field 

2 

dt{pO) = --pdaiipapa (28) 

We now make the Ansatz 

a a p(9) = p9a a(3 (29) 
so that we get from (^) for the time dependence of the temperature 

0{t) = 6 Q cxp (^^d a u l 3CF f 3 a t S j (30) 

which reads for the velocity field of eqn. ( p3| ) 

6(t) = 6 exp(-p& X yt \ . (31) 

So the differential equations for a become algebraic equations 

Oxx + r ( -^ja xy ) a xx = r ( ~ - 2 ) jcr xy (32) 



Gyy +T 



2 ' y j ■**«*■ i ^ 

2 



jVxy Oyy = T~ja x y (33) 
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n ,n - r I ~^rjo- xy I cr^y + r-ycr^y = -rp9j (34) 
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Thus we arrive at a cubic equation for a xy 



- —<?xy - 1 ) cr xy = -rpdj (35) 



2 T7 

3 ~pl 

The solution can be written in analytic terms for r = tq j p as 

9{t) = 9 exp(t\) (36) 
/ 1 2t 2 7 2 \ 

°-" = MiTTa + oTTaf) (37) 

*vv = P ®!^ ( 38 ) 

(4°) 
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with 



4 . .o/l 



A = — sinh 2 I - cosh" 1 (1 + 9t 2 Y) ) (41) 



which is also the exact solution for the BGK approximation of the Boltzmann 
equation first described by Zwanzig || where he derived /„ (f ) from eqn. (f| for 
a simple shear flow. He then used this solution in the equivalent of eqn. (@) to 
determine the stress. We have neglected transient terms which are decaying ex- 
ponentially. We recover the exact result because our approximation of eqn.([j"7|) 
become exact for spatially constant temperature. 

There are three invariants of the stress tensor which are defined as a xy = —rfy 
related to the viscosity 77, a xx — a yy = — ^ij 2 related to the first normal stress 
coefficient and a yy — a zz = — V E , 27 2 related to the second normal stress 
coefficient ^2- For our constitutive equation we get 

In this simple situation we get a result that is equivalent to that of an up- 
per convected Maxwell model but with the opposite sign for the first normal 
stress difference. Viscoelastic polymeric systems for which the upper convected 
Maxwell model was devised seem to always have a positive first normal stress 
coefficient v&i. In the next section we will give some intuitive understanding for 
this fundamental difference in the viscoelasticity. 

Although for a simple shear flow no stationary solution exists due to viscous 
heating a stationary state can be achieved with the help of a thermostat. The 
simplest way of implementing a thermostat is simply to assume that the partial 
time derivatives and the spatial derivatives of the temperature in eqns (^5|)-(^8|) 
vanishes (|, This is equivalent to a situation where a thermostat insures 
a constant temperature, but does not otherwise influence the dynamics. This 
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Figure 1: The exact solution for the shear-rate dependent viscosity 77(7) and 
first normal stress difference against the shear rate 7 (full line) and the 

isothermal solution (broken line). See text for details. 



thermostat will change the energy equation but not the constitutive equation. 
In this case we get 

V= : t T f 2 . 2 *i = -2t»j * 2 = (43) 
1 + |t^7 2 

and we see that the shear-thinning persists although it now has a different form. 
The first normal stress difference now only depends on the shear rate through the 
viscosity. The viscosity and the first normal stress difference for the isothermal 
and the heating case are compared in Figure [|. Other thermostats lead to 
different results. Lee and Dufty [J considered a thermostat that is frequently 
used for Molecular Dynamics simulations. This thermostat rescales the local 
velocities relative to the average velocities such that the temperature remains 
constant. In this case the constitutive equation is also changed but in a way 
that recovers the same solution [l2] as in the case with heating, except that the 
temperature remains constant. 



3.2 Extensional flow 

We will now consider a second family of flow-fields for which, to my knowledge, 
the BGK equations had not previously been solved. We will, again, assume a 
spatially constant temperature so that the results will be exact. And extensional 
flow is a shear-free flow defined by the velocity profile 

/ -(l + b)x/2 \ 
u = e I -(l-6)y/2 j (44) 

where < b < I and e is the elongation rate. Several special shear-free flows 
are obtained for particular choices of the parameter b: 
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Elongational flow: (6 — 0, e > 0) 

Biaxial stretching flow: (6 — 0, e < 0) 

Planar elongational flow: (6 = 1) 

For the time-dependence of the temperature we get from eqn. (|3^) for the 
extensional flow 

9{t) = 9 exp + b)d xx - (1 - 6)<r ra + 2a zz )t \ . (45) 

The flow is a linear flow and if we assume a constant initial temperature 
the solutions of the constitutive equation are again exact. For the constitutive 
equation for the stress (|l^), using that <r xx + <r yy + a zz — 0, we get 



(-L - (1 + b))a xx + Ul + ~)a xx + (1 - ~)a yil ) ( i + <r r i = I + b ( 4(,) 

(A - (1 - 6))5- ra + ((1 + ^)& xx + (1 - (1 + = 1-6 (47) 

These equations have analytic solutions one of which is physical. These results, 
however, are too lengthy to be presented here. For the simpler case of b = 
the x and y directions are symmetric and so we can impose a xx = a yy and the 
equations can be simplified to 

(^t - \)o xx + 2a xx (1 +a xx ) = 1 (48) 
re 

This equation has two solutions but one of these solutions is unphysical as the 
stress is larger that pQ which would require a negative probability density /„ in 
(||). The physical solution is 



pOl + re- v / l + 2re + 9r 2 e 2 

Oxx = ~~ A : (49) 

4 re 

and ~a zz /2 — a yy — a xx . 

It is customary to define two viscosity functions to describe the rheological 
behaviour of a fluid in extensional flow as 

<Jxx - <J ZZ = fji(e,b)e (50) 
o- xx -<J yy = mie^fye (51) 

and we show the values for different b in Figure |^. For the elongational flow 
and the biaxial stretching flow (b = 0) defined above the x and y-directions are 
equivalent and therefore fj 2 — 0. For the elongational flow (b — 0, e > 0) we find 
only shear thinning (in polymeric flows a shear-thickening is observed) but for 
the biaxial stretching flow gases show a shear thickening. 

For the planar elongational flow 6=1 the x and z directions become sym- 
metric if you invert the extension rate (a xx (e, 6=1) = a zz (— e, 6 = 1)) so that 
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TC 

Figure 2: The exact solution for the extension-rate dependent elongational vis- 
cosities 771 (full lines) and ^(dashed lines) for different parameters b (see text 
for details). 



the viscosity 771 becomes symmetric, but the viscosity 772 becomes asymmetric. 
Even though the planar elongational flow is a two-dimensional flow the stress 
will be different for a true two-dimensional system and this three dimensional 
case. 

If we consider a system with a simple-minded thermostat again that simply 
assumes a constant temperature we get unphysical divergences for the stress 

fl . 3(l + b) + r6(-3 + 6& + 6 2 ) 
"~ = ^ e -3 + (3 + * )T 2*2 (52) 

Q . 3(l-6) + T6(-3-6b + 6 2 ) 

= Tp9e -3 + (3 + y)r^ (53) 

o n- -3 + re(3-fo 2 ) 

but the quadratic behaviour for small e still agrees with the non-isothermal case. 
These divergences indicate that the simple-minded thermostat is unphysical. 

It is easy to see why this unphysical behavour can occur. Because the ther- 
mostat will regulate the temperature without influencing the stress we can have 
a second moment 

n_ = //„(*-«.)" (55) 

that becomes negative which in terms reqires negative contributions to the prob- 
ability distribution. This is what makes the thermostat unphysical. 



4 Intuition for the viscoelastic effects in a gas 

In Figures || and ^ we sketched the origins of the non-isotropic stress for both the 
case of a gas and the case of a polymeric liquid in a simple shear flow. In a gas 
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(a) 



(b) 



Figure 3: Illustration of the origin of a non-diagonal pressure II and the corre- 
sponding stress tensor a in the case of a gas. (a) Origin of pressure tensor (b) 
the traceless stress tensor (vectors pointing inwards represent a negative stress) . 
Note that the orientation of the stress tensor we have 8 n > 45o. 




Figure 4: Sketch of the origin and orientation of the non-diagonal stress tensor 
for a polymeric liquid in a simple shear flow. The stress looks qualitatively 
similar to the gas case of Figure || but note that here 6 P < 45°. Also there is no 
requirement for the polymeric stress to be traceless. 
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the mean velocity of the particles is given by the flow profile of equation (g3|) . 
Particles that are convected in the direction of —y carry more x-momentum 
with them than the average at this position and particles that are convected 
in the direction of +y carry less x-momentum with them that the average that 
the new position. Therefore the momentum distribution will be non-isotropic 
as indicated in Figure ^(a). If we subtract the isotropic part of the momentum 
distribution we get the stress tensor a of equation (|8|) shown in Figure ||(b). Note 
that the angle 9 g defined in this Figure is always smaller that 45°. Therefore the 
non-equilibrium stress tends to reduce the force on the walls which is equivalent 
to a negative first normal stress coefficient. 

In simple shear flow the viscosity is a measure for the transport of x- momentum 
in the y-direction. Because there are now fewer particles streaming in the y- 
directions this also means that the viscosity is reduced. This is the intuitive 
reason for the shear thinning in a gas. 

In the case of a polymeric liquid the origin of the non-isotropic stress lies 
in the stretching of the macromolecules as sketched in Figure |J. In equilibrium 
without flow the macromolecules would curl up into a spherical shape but the 
flow tends to deform the shape to a more ellipsoidal from. The result is a 
restoring force that would restore the molecule to a spherical shape. The first 
deformation of the molecule occurs at 9 p = 45° and increases from there on. The 
effect is that the stress distribution sketched in Figure [I] increases the pressure 
in the y-direction relative to the pressure in the x-direction which is equivalent 
to a positive first normal stress coefficient. 

The deformation of the coils also means that there is less transport of x- 
momentum in the y-direction which is the reason for shear-thinning in the poly- 
meric liquid. 

5 Conclusions 

We have derived a viscoelastic constitutive equation for a gas described by the 
BGK approximation of the Boltzmann equation and shown that this constitutive 
equation differs substantially from those that describe the viscoelastic properties 
of polymeric materials. For simple shear flow the exact results for the BGK 
approximation of the Boltzmann equation are recovered. We explained the 
different sign of the first normal stress difference in gases and polymeric liquids 
be examining the qualitative different origins of the non-Newtonian stress. 
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